;************************************************
; NCL Graphics: 
;************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" 
load "$NCARG_LIB/ncarg/nclscripts/csm/contributed.ncl"

;************************************************

begin

  load "load_setup.ncl"

  vname  = "rh"
  vunits = "K"
  vscale = 1. ;;e-3

  pa = basepa+ca+"/"
  pb = basepb+cb+"/"
  fa = ca+ma
  fb = cb+mb
  caseM  = (/pa+fa,pb+fb/)
  case  = (/caseM(0),caseM(1)/)
  Ncase = dimsizes(case)
  plotname = "figure_"+vname+"_zonal_"+shortname

  print(" ")
  print(" ")
  print("..................................................................")
  print(" ")
  print(" ")
  print(" varname   : "+plotname)
  print(" ctrl      : "+case(0))
  print(" test      : "+case(1))
  print(" plot name : "+plotname)
  print(" ")
  print(" ")

;************************************************
; set parameter
;************************************************

  caseLP  = (/pa+fa/) 
  caseBN  = (/pb+fb/) 
  Ncase   = 2 
  
  labelsLP=(/"MOZAIC","AIRS","CNT","CNT max 50%"/) 
  Nlabels = 4 

  
  title = "" 
  
  Nbin=100
    
  number=(/"0","1","2","3","4","5","6","7","8","9"/)
  
;************************************************
; read in obs variables
;************************************************

  obs1=asciiread("mozaicRH_fmt_RHpdf.txt",(/Nbin,2/),"float")
  obs2=asciiread("airsRH_RHpdf.txt",(/Nbin,2/),"float")
  ;printVarSummary(obs)
  
;************************************************
; read in obs variables
;************************************************

   filename=caseLP(0)
   fi=addfile(filename,"r")
   lev   = fi->lev
   lat   = fi->lat
   lon   = fi->lon 
   nlev  = dimsizes(lev)                           
   nlat  = dimsizes(lat)
   nlon  = dimsizes(lon)
   ;print("nlev "+sprinti("%7.2i",nlev))
   ;print("nlat "+sprinti("%7.2i",nlat))
   ;print("nlon "+sprinti("%7.2i",nlon))
   
   
   fre = new((/Ncase,Nbin/), float)
   do n=0,1
    if (n.eq.0) then
       case=caseLP
    else
       case=caseBN
    end if
    do f=0,Ncase-1
     filename=case(f)
     ;print ("open file :: "+ filename)
     ;print ("read form  "+filename)
     fi=addfile(filename,"r") 
     do n1=0,9
     do n2=0,9
     bin=n1*10+n2
     varnames="RHi"+number(n1)+number(n2)
     Vartema = dim_avg_Wrap(fi->$varnames$(time|:,lev|:, lat|:, lon|:))
     Vartemp = Vartema(0,:,:) 
     ;print(" "+varnames)
     count=0.
     fre(f,bin)=0.
     do nz=0,nlev-1
     if ((lev(nz).ge.150).and.(lev(nz).le.300)) then
        do ny=0,nlat-1
        if ((lat(ny).ge.-30).and.(lat(ny).le.60)) then
           fre(f,bin)=fre(f,bin)+Vartemp(nz,ny)
           count=count+1.
        end if
        end do
     end if
     end do
     if (count.ge.1) then
        fre(f,bin)=fre(f,bin)/count
     else
        fre(f,bin)= -9999.99
     end if
     end do
     end do
    end do
    freLP=fre
   end do
  
;************************************************
; variables for plot
;************************************************

  YLP= new((/Nlabels,Nbin*2/), float)
  X  = new(          Nbin*2  , float)
  X(0)        = 0
  X(Nbin*2-1) = 200
  do i=1,Nbin-1
     X(i*2)  =2*i
     X(i*2-1)=2*i
  end do
  do i=0,Nbin-1
     if (obs1(i,0).lt.1e-5) then 
        obs1(i,0)=1e-5
     end if
     YLP(0,i*2  )=obs1(i,0)
     YLP(0,i*2+1)=obs1(i,0)
     if (obs2(i,0).lt.1e-5) then 
        obs2(i,0)=1e-5
     end if
     YLP(1,i*2  )=obs2(i,0)
     YLP(1,i*2+1)=obs2(i,0)
     do f=0,Ncase-1
        if (freLP(f,i).lt.1e-5) then 
           freLP(f,i)=1e-5
        end if
        YLP(f+Nlabels-Ncase,i*2  )=freLP(f,i)
        YLP(f+Nlabels-Ncase,i*2+1)=freLP(f,i)
     end do
  end do
  
;************************************************
; set default res
; rea2 for gsn_csm_xy
;************************************************
  wks   = gsn_open_wks (plotfmt,plotname)    ; open workstation
  gsn_define_colormap(wks,"gui_default")          ; choose colormap
  res                      = True                ; plot mods desired
  res@gsnFrame             = False               ; don't advance page
  res@gsnDraw              = False               ; don't advance frame
  res@vpXF                  = 0.20     ; x location left
  res@vpYF                  = 0.80     ; y location  up
  res@vpWidthF              = 0.65     ; width
  res@vpHeightF             = 0.40     ; height
  res@gsnStringFontHeightF  = 0.02
  res@tiMainString          = "  "    ; title  
  res@gsnLeftString = "RH PDF"
  ;res@gsnRightString= "30S-60N 200-500hpa"
  res@gsnRightString= "30S-60N 150-300hpa"
  res@tiXAxisString = "RH(%)"
  res@tiYAxisString = "Relative frequency"
  res@xyYStyle      = "Log"            ;Log, Linear, Irregular 
  res@xyComputeXMin = True 
  res@xyLineColors      = (/"red","RoyalBlue","black", "LimeGreen", "DarkSalmon"/)  ; line colors  
  res@xyLineThicknesses = (/1.5,1.5,2.0,2.5,2.5/)        ; line thicknesses
  res@xyDashPatterns    = (/0  ,0  ,0  ,0  ,0  /)        ; line patterns
  res@trYMinF           = 1e-5             ; Y axis min 
  res@trYMaxF           = 1.0              ; Y axis max
  res@trXMinF           = 0                  ; X axis min 
  res@trXMaxF           = 200                ; X axis max
  res@tmXBMode = "Explicit"
  res@tmXBValues = (/0,50,100,150,200/)
  res@tmXBLabels = (/"0","50","100","150","200"/) 
  res_text                    = True                  ; text mods desired
  res_text@txFontHeightF      = 0.021                 ; change text size
  res_text@txJust             = "CenterLeft"          ; text justification
  res_lines                   = True                  ; polyline mods desired
  res_lines@gsLineDashPattern = 2.                    ; solid line
  res_lines@gsLineThicknessF  = 0.7                   ; line thicker
  resP                 = True          
  resP@gsnPanelYWhiteSpacePercent = 0  ; draw panel with white space added
  resP@gsnPanelXWhiteSpacePercent = 2  ; draw panel with white space added
  resP@txString            = title 
  
;************************************************
; create plot
;************************************************   

  plot = gsn_csm_xy (wks,X,YLP,res)
  
  xx = (/100         ,100 /)
  yy = (/res@trYMinF,res@trYMaxF/)
  addline0 = gsn_add_polyline(wks,plot,xx,yy,res_lines)             ; add polyline

  xxx=30
  yyy=8e-3
  do v=0,Nlabels-1
     yyy=yyy/exp(0.9)
     res_text@txFontColor = res@xyLineColors(v)
     addtext0=gsn_add_text(wks,plot,labelsLP(v),xxx,yyy,res_text)
  end do

  gsn_panel(wks,plot,(/1,1/),resP)
 
  ;print("ps file :  "+title)
  
end



